NASA/CR-2002-211655 
ICASE  Report  No.  2002-15 


Covolume-based  Intergrid  Transfer  Operator  in  Pj 
Nonconforming  Multigrid  Method 


Kab  Seok  Kang 
ICASE,  Hampton,  Virginia 


May  2002 


Report  Documentation  Page 


Report  Date 

Report  Type 

Dates  Covered  (from...  to) 

00MAY2002 

N/A 

- 

Title  and  Subtitle 

Contract  Number 

Covolume-based  mtergrid  Transrer  Operator  in  PI 
Nonconforming  Multigrid  Method 

Grant  Number 

Program  Element  Number 

Author(s) 

Project  Number 

Kab  Seok  Kang 

Task  Number 

Work  Unit  Number 

Performing  Organization  Name(s)  and  Address(es) 

ICASE  Mail  Stop  132C  NASA  Langley  Research  Center 
Hampton,  VA  23681-2199 

Performing  Organization  Report  Number 

Sponsoring/Monitoring  Agency  Name(s)  and 

Sponsor/Monitor’s  Acronym(s) 

Address(es) 

Sponsor/Monitor’s  Report  Number(s) 

Distribution/Availability  Statement 

Approved  for  public  release,  distribution  unlimited 

Supplementary  Notes 

ICASE  Report  No.  2002-15 

Abstract 

see  report 

Subject  Terms 

Report  Classification 

unclassified 

Classification  of  this  page 

unclassified 

Classification  of  Abstract 

unclassified 

Limitation  of  Abstract 

SAR 

Number  of  Pages 

19 

The  NASA  STI  Program  Office  ...  in  Profile 


Since  its  founding,  NASA  has  been  dedicated 
to  the  advancement  of  aeronautics  and  space 
science.  The  NASA  Scientific  and  Technical 
Information  (STI)  Program  Office  plays  a  key 
part  in  helping  NASA  maintain  this 
important  role. 

The  NASA  STI  Program  Office  is  operated  by 
Langley  Research  Center,  the  lead  center  for 
NASA’s  scientific  and  technical  information. 
The  NASA  STI  Program  Office  provides 
access  to  the  NASA  STI  Database,  the 
largest  collection  of  aeronautical  and  space 
science  STI  in  the  world.  The  Program  Office 
is  also  NASA’s  institutional  mechanism  for 
disseminating  the  results  of  its  research  and 
development  activities.  These  results  are 
published  by  NASA  in  the  NASA  STI  Report 
Series,  which  includes  the  following  report 
types: 

•  TECHNICAL  PUBLICATION.  Reports  of 
completed  research  or  a  major  significant 
phase  of  research  that  present  the  results 
of  NASA  programs  and  include  extensive 
data  or  theoretical  analysis.  Includes 
compilations  of  significant  scientific  and 
technical  data  and  information  deemed 

to  be  of  continuing  reference  value.  NASA’s 
counterpart  of  peer-reviewed  formal 
professional  papers,  but  having  less 
stringent  limitations  on  manuscript 
length  and  extent  of  graphic 
presentations. 

•  TECHNICAL  MEMORANDUM. 

Scientific  and  technical  findings  that  are 
preliminary  or  of  specialized  interest, 
e.g.,  quick  release  reports,  working 
papers,  and  bibliographies  that  contain 
minimal  annotation.  Does  not  contain 
extensive  analysis. 

•  CONTRACTOR  REPORT.  Scientific  and 
technical  findings  by  NASA-sponsored 
contractors  and  grantees. 


•  CONEERENCEPUBEICATIONS. 
Collected  papers  from  scientific  and 
technical  conferences,  symposia, 
seminars,  or  other  meetings  sponsored  or 
cosponsored  by  NASA. 

•  SPECIAE  PUBEICATION.  Scientific, 
technical,  or  historical  information  from 
NASA  programs,  projects,  and  missions, 
often  concerned  with  subjects  having 
substantial  public  interest. 

•  TECHNICAE  TRANSEATION.  English- 
language  translations  of  foreign  scientific 
and  technical  material  pertinent  to 
NASA’s  mission. 

Specialized  services  that  complement  the 
STI  Program  Office’s  diverse  offerings  include 
creating  custom  thesauri,  building  customized 
data  bases,  organizing  and  publishing 
research  results  .  .  .  even  providing  videos. 

Eor  more  information  about  the  NASA  STI 
Program  Office,  see  the  following: 

•  Access  the  NASA  STI  Program  Home 
Page  at  http://www.sti.nasa.gov 

•  Email  your  question  via  the  Internet  to 
help@sti.nasa.gov 

•  Eax  your  question  to  the  NASA  STI 
Help  Desk  at  (301)  621-0134 

•  Telephone  the  NASA  STI  Help  Desk  at 
(301)  621-0390 

•  Write  to: 

NASA  STI  Help  Desk 

NASA  Center  for  AeroSpace  Information 

7121  Standard  Drive 

Hanover,  MD  21076-1320 


NASA/CR-2002-211655 
ICASE  Report  No.  2002-15 


Covolume-based  Intergrid  Transfer  Operator  in  Pj 
Nonconforming  Multigrid  Method 


Kab  Seok  Kang 
ICASE,  Hampton,  Virginia 


ICASE 

NASA  Langley  Research  Center 
Hampton,  Virginia 

Operated  by  Universities  Space  Research  Association 


Prepared  for  Langley  Research  Center 
under  Contract  NAS  1-97046 


May  2002 


Available  from  the  following: 

NASA  Center  for  AeroSpace  Information  (CASI) 
7121  Standard  Drive 
Hanover,  MD  21076-1320 
(301)  621-0390 


National  Technical  Information  Service  (NTIS) 
5285  Port  Royal  Road 
Springfield,  VA  22161-2171 
(703)  487-4650 


COVOLUME-BASED  INTERGRID  TRANSFER  OPERATOR  IN  Pi  NONCONFORMING 

MULTIGRID  METHOD 


KAB  SEOK  KANG* 

Abstract.  In  this  paper,  we  introduce  an  intergrid  transfer  operator  which  is  based  on  the  covolume 
of  nodes  in  a  Pi  nonconforming  multigrid  method  and  study  the  convergence  behavior  of  the  multigrid 
method  with  this  intergrid  transfer  operator.  This  intergrid  transfer  operator  needs  fewer  computations  and 
neighborhood  node  values  than  previous  operators,  which  is  a  good  property  for  parallelization.  The  Pi 
nonconforming  multigrid  method  with  this  intergrid  transfer  operator  is  suitable  for  solving  problems  with 
Robin  boundary  conditions  and  nonlinear  problems  with  bound  constraints  on  solutions. 

Key  words,  multigrid  method,  covolume  method,  nonconforming  finite  elements,  elliptic  equations 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Multigrid  methods  are  well  known  as  efficient  solution  techniques  for  many  prob¬ 
lems  including  elliptic  and  hyperbolic  partial  differential  equations,  nonlinear  problem,  and  even  systems 
of  algebraic  equations  that  are  not  derived  from  the  spatial  discretization  of  a  partial  differential  equation 
([5,  10,  14,  17]).  Because  nonconforming  finite  element  or  covolume  methods  have  proven  flexible  and  effec¬ 
tive  on  incompressible  fluid  flow  problems  and  biharmonic  and  plate  problems  ([9,  8,  11])  many  researchers 
have  been  interested  in  studying  multigrid  methods  for  nonconforming  finite  elements  or  covolume  methods 
([6,  7,  2,  4,  1,  12]). 

In  nonconforming  multigrid  methods,  the  intergrid  transfer  operators  have  important  roles  in  conver¬ 
gence.  In  this  paper,  we  consider  a  covolume-based  intergrid  transfer  operator.  This  intergrid  transfer 
operator  needs  less  computation  and  neighborhood  node  information  than  previously  proposed  intergrid 
transfer  operators.  The  Pi  nonconforming  multigrid  method  with  previous  intergrid  transfer  operators  is 
less  suitable  for  solving  nonlinear  problems  that  have  bounds  on  solutions  and  has  poor  error  reduction 
for  problems  with  Robin  boundary  conditions.  However  the  Pi  nonconforming  multigrid  method  with  this 
intergrid  transfer  operator  is  very  suitable  for  such  problems. 

Many  authors  have  shown  that  nonconforming  TT-cycle  multigrid  methods  converge  and  nonconforming 
variable  U-cycle  multigrid  preconditioners  have  a  uniform  condition  number  as  preconditioners.  In  this 
paper,  we  investigate  the  convergence  behavior  of  hh-cycle  multigrid  methods  and  the  condition  number  as 
preconditioners  of  U-cycle  multigrid  methods  with  covolume-based  intergrid  transfer  operators,  by  means  of 
numerical  experiments. 

This  paper  is  organized  as  follows.  In  section  2,  we  summarize  some  results  of  Pi  nonconforming  finite 
element  and  covolume  methods.  In  section  3,  we  introduce  the  covolume-based  intergrid  transfer  operator 
and  recount  the  abstract  theory  developed  by  Bramble  et  al.  for  nonnested  multigrid  methods.  In  section 
4,  we  report  the  results  of  numerical  experiments  justifying  the  convergence  theory  presented  in  section  3 
and  applied  to  a  Radiation  Transport  problem  which  is  a  system  of  coupled  nonlinear  partial  differential 
equations  with  discontinuous  diffusion  coefficients. 

*ICASE,  Mail  Stop  132C,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199  email:  kks002@icase.edu.  This 
research  was  partially  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract  No.  NASl- 
97046  while  the  author  was  in  residence  at  ICASE,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199.  This  work  was 
also  partially  supported  by  postdoctoral  fellowships  program  from  Korea  Science  &  Engineering  Foundation  (KOSEF). 
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2.  Model  problem  and  its  discretizations.  We  consider  the  second-order  elliptic  problem  with 
Robin  boundary  conditions 


—V  •  AVu  =  /,  in  fl, 

l3u  +  A—  =  o,  on  dfl, 
on 


(2.2.1) 


where  is  a  bounded  polygonal  domain  with  boundary  dfl,  /  €  L‘^{fl),  g  G  L‘^{dft),  and  the  symmetric 

coefficient  matrix  A  €  satisfies 


aox\  <  X^Mx,  y)x  <  {x,  y)  &  fl,  X  G  IR^ 


(2.2.2) 


with  fixed  constants  a(y,ai  >  0.  It  is  well  known  that  if  (3  is  not  equal  to  zero  on  some  set  of  boundary 
which  is  not  of  measure  zero,  then  equation  (2.2.1)  has  a  unique  solution. 

Let  H^{ft)  and  L‘^{ft)  be  the  usual  Sobolev  spaces  with  norm  and  seminorm 


where  s  is  a  nonnegative  integer.  Let  (•,•)  denote  the  L‘^{ft)  inner  product.  As  usual,  the  L‘^{ft)  norm  is 
indicated  by  ||  •  ||o. 

The  variational  form  of  (2.2.1)  can  be  written  as  follows:  Find  u  G  H^{ft)  such  that 

a{u,v)  =  {f,v)  +  {g,v)dQ,  \/v  e  H^(fl),  (2.2.3) 

where  a{v,w)  =  (AVVjVw)  +  {f3v,w)QQ,  for  all  v,w,  G  H^{ft). 

Let  /lo  and  Thq  =  To  be  given,  where  To  is  a  partition  of  ft  into  triangles  and  /lo  is  the  maximum  diameter 
of  the  elements  of  To-  For  each  integer  1  <  k  <  J ,  lei  hu  =  2“*/io  and  the  sequence  of  triangulations 
Th^  =  Ti;  be  constructed  by  the  nested-mesh  subdivision  method,  i.e.,  let  Tk  be  constructed  by  connecting 
the  midpoints  of  the  edges  of  the  triangles  in  Ti;-i,  and  lei  Th  =  Tj  be  the  finest  grid. 

Define  the  Pi -nonconforming  finite  element  spaces 

14  =  {u  G  L‘^{ft)  :  u lit  is  linear  for  all  IL  G  T*, 

V  is  continuous  at  the  midpoints  of  interior  edges}. 


In  Pi -nonconforming  finite  elements,  the  node  points  are  the  midpoints  of  the  edges.  Obviously,  this  results 
in 


14  sT  14  14  =  14. 

Define  the  bilinear  forms  over  the  spaces  14,  for  each  k,  as  follows: 

ak{v,w)  =  {AVv,Vw)k  +  {v,w)8n,  \/v,wGVk, 

Ke% 

where  (•,  ■)k  is  the  inner  product.  Moreover,  define  the  discrete  energy  norms  as  follows: 

\\v\\i,k  =  ak{v,vY^‘^ ,  VuG14- 
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Fig.  1.  Primal  and  Dual  element 

Then  the  nonconforming  finite  element  discretization  equation  of  (2.2.3)  can  be  written  as:  Find  Uh  G  Vh 
such  that 


ah{uh,Vh)  =  {f,Vh)  +  {g,Vh)dii,  Vh&Vh-  (2.2.4) 

From  a  previous  study  of  Pi  nonconforming  finite  element  methods  [9] ,  we  have 

\\u  -  UhWo  +  h\\u  -  Uh\\i,k  <  Chl\u\2  (2.2.5) 

if  u  G 

To  define  a  covolume  of  nodes,  we  construct  the  dual  partitions  .  Divide  each  triangle  of  the  primal 
partition  into  three  sub-triangles  by  connecting  two  vertices  and  the  barycenter  of  a  primal  element  as  in 
Figure  1(a).  As  in  Figure  1(b),  the  dual  element  based  at  the  node  Pi  (covolume  of  Pi)  is  made  up  of  the  two 
triangles  AAiCiA^  and  AAiC^A^.  We  do  the  obvious  modification  at  a  boundary  node.  Carrying  out  the 
construction  for  every  node  in  the  primal  partition,  we  obtain  a  dual  partition  for  the  domain.  We  denote 
the  covolume  of  node  P  as  Kp  and  the  dual  partition  as  =  UKp.  Define  the  associated  test  function 
spaces  Yk  as  the  space  of  piecewise  constant  functions: 

Yk  =  {z&  L^n)  :  z\k^  is  a  constant  vector}. 

Obviously,  we  have 

Yo<tYi<t---<tYj=Yh. 


Define  an  operator  from  the  spaces  Vk  xYk,  for  each  k,  as  follows. 


a*{uh,Vh) 


Nu 

E 


~  JdKi,,\dKi,,ndn  un 

.OUh 


fduhVhda 


9Q 


J2vh{Pi) 


i^l 


dK^,\dK^,ndQ 


dn 


I  . 


da  jduhVhda^ 


where  is  the  outer  normal  derivative  of  Uh- 

The  covolume  discretization  equation  of  (2.2.1)  can  be  written  as:  Find  u\  G  Vh  such  that 


aliul^Vh)  =  {f,Vh)  +  {g,Vh)dii,  Vuft  G  Yh. 


(2.2.6) 


(2.2.7) 
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Q  :  node  of  coarse  grid 
•  :  node  of  fine  grid 

Fig.  2.  Covolume  of  coarse  and  fine  grid 


We  introduce  one-to-one  transfer  operator  7^;  from  V*  54  by 

ikVhix)  =  ^vh{Pj)xj{x),  yx  G  n, 

i=i 

where  Xj  is  a  characteristic  function  associated  with  the  dual  element  Kfi. ,  j  =  1, . . .  ,Nk. 
If  we  assume  that  A  is  piecewise  constant,  then  we  have  the  following  relation 

ak{uh,Vh)  =  al{uh,XkVh),  yuh,Vh  G  Vk 


from  [8]. 

We  have  the  following  convergence  estimate  from  [8]  and  [11]: 


w  “  Wftlli.i:  <  C/idlwIb  +  1),  if  u  G 
\W  ~  Wftllo  <  +  1),  if  u  G 


(2.2.8) 


(2.2.9) 


(2.2.10) 


In  the  case  where  A  is  piecewise  constant,  the  finite  element  method  and  the  covolume  method  differ 
only  in  the  source  term  (right-hand  side  of  equation).  Therefore,  in  the  next  section,  we  consider  only  the 
discretization  of  the  linear  problem  induced  by  the  finite  element  method. 

3.  Multigrid  method  and  covolume-based  intergrid  transfer  operator.  In  this  section,  we 
introduce  a  covolume-based  intergrid  transfer  operator  and  abstract  theories  developed  by  Bramble  et  al. 
for  non-nested  multigrid  algorithms. 

The  covolume-based  intergrid  transfer  operators  It  ■  Vu-i  V*  for  A:  =  1, . . .  ,  J  are  defined  by 

Nu-i 

Ikv{P)  =  ^  v{Pi)\K*p^  n  K*p\/\K*p\,  (3.3.1) 

where  P  is  a  node  point  of  Tk  and  |A|  is  the  area  of  A. 
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Remark  3.1.  We  consider  Figure  2  as  an  example.  Then 

ihvm)  =  and  (hv)(Pi)  =  v(Pi). 

Compared  to  other  intergrid  transfer  operators  in  [2](ll)  and  [ll](ll^), 

=  '^{Pi)  +  ■^{^{Pii)  +  v(Piii)  -  v(Piv)  -  v(Pv)} 

(4"^)m)  =  \v{Pi)  +  ^{v{Pu+v{Pui)}, 

this  operator  needs  less  computation  to  compute  values  of  nodes  that  are  located  in  edges  of  the  coarse  grid. 

Remark  3.2.  The  intergrid  transfer  operator  proposed  in  [2]  cannot  be  employed  in  some  nonlinear 
problems  that  have  a  restricted  range  of  solution  values  and  rapidly  changing  fields,  such  as  the  Marshak 
wave  problem  in  radiation  transport.  If  values  of  the  solution  are  everywhere  positive  and  v(Pi),  v(Pii)  and 
v{Piii)  are  small  and  v{Piv)  and  v{Pv)  are  big,  then  {llv){Pi)  is  negative. 

Let  Ak  :Vk  ^  14  be  the  discretization  operator  on  level  k  given  by 

(AkV,w)k  =  ak(v,w),  \/v,weVk. 

The  operator  A^  is  clearly  symmetric  (in  both  the  ak{-,-)  and  (•,  •)  inner  products)  and  positive  definite. 
Also,  we  define  the  operators  P^_i  ■  14  ^  14-i  and  Pk-i  ■  14  ^  14-i  by 

(IkV,w)  =  (u,P°_iw),  Vu  G  14_i,Vw  G  14, 


and 


ak-i{Pk-iw,v)  =  OkiwJkv),  \fv  G  14_i,Vw  G  14. 

It  is  easy  to  see  that  IkPk-i  is  a  symmetric  operator  with  respect  to  the  Ok  form.  Note  that  neither 
nor  i4_i  is  a  projection  in  the  nonconforming  case. 

Remark  3.3.  Because,  for  a  coarse  triangulation  node  P,  Kp,  n  Kp  =  0  except  for  six  covolumes  Kp. 
maximally  in  a  fine  triangulation,  the  number  of  fine  triangulation  nodal  values  which  need  to  calculate  Pk-i 
is  less  than  six  and  this  is  less  than  the  number  of  required  fine  triangulation  nodal  values  in  the  intergrid 
transfer  operator  previously  proposed  in  [2],  [12].  As  an  example,  we  consider  Figure  2.  The  transfer  operator 
which  was  proposed  in  [2]  requires  the  values  of  Pi,. . .  ,  P14;  and  the  transfer  operator  which  was  proposed  in 
[11]  requires  the  values  of  Pi,  ■  ■  ■  ,  T’lo  to  compute  the  value  of  P.  But  the  transfer  operator  (3.3.1)  needs  only 
Pi,...  ,Pq.  Because,  in  the  computation  of  Ik  and  Pk-i,  we  need  nodal  values  to  calculate  multiplication 
by  the  matrix  Ak,  this  property  is  good  for  parallelization  when  the  domain  is  partitioned  by  cutting  edges 
between  vertices  and  to  preserve  covolumes,  so  each  vertex  is  uniquely  owned. 

Finally,  let  :  14  14  for  A:  =  1, . . .  ,  J  be  the  linear  smoothing  operators,  let  denote  the  adjoint 

of  Rk  with  respect  to  the  (•,  •)  inner  product,  and  define 

^(i)  _  \  Rk,  ^  odd, 

[R[  ,  I  even. 

Following  [3],  the  multigrid  operator  Bk  -Vk  ^  14  is  defined  recursively  as  follows. 

Multigrid  Algorithm  3.1.  Let  1  <  A:  <  J  and  p  be  a  positive  integer.  Set  Bq  =  A((^ .  Assume  that 
Bk-i  has  been  defined  and  define  Bkg  for  gi  G  14  by 
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(1)  Set  =  0  and  =  0. 

(2)  Define  for  f  =  1, . . .  ,  m(k)  by 

x‘=x‘-^+Rt"^'^^\9-A,x‘-^). 

(3)  Define  =  a;™(*)  -|-  where  g*  for  i  =  1, . . .  ,p  is  defined  by 

g'  =  g'-'  +  Bk-i[Pt,{g  -  -  Ak-W-^]. 

(4)  Define  for  I  =  m(k)  +  1, . . .  ,  2m{k)  by 

(5)  Set  Bug  = 

In  Multigrid  algorithm  2.1,  m{k)  gives  the  number  of  pre-  and  post-smoothing  iterations  and  can  vary 
as  a  function  of  k.  If  p  =  1,  we  have  a  F-cycle  multigrid  algorithm.  If  p  =  2,  we  have  a  W-cycle  multigrid 
algorithm.  Other  versions  of  multigrid  algorithms  without  pre-  or  post-smoothing  iterations  can  be  analyzed 
similarly.  A  variable  F-cycle  multigrid  algorithm  is  that  for  which  the  number  of  smoothing  m{k)  increases 
exponentially  as  k  decreases  (i.e.,  p  =  1  and  m{k)  =  2'^“*). 

Based  on  the  methodology  developed  in  [3] ,  two  other  very  important  ingredients  in  convergence  analysis 
of  non-nested  multigrid  algorithms  are  the  following  assumptions. 

A.l:  (Regularity  and  approximation  assumption)  For  some  a  €  (0,1]  there  exists  Ca,  independent  of 
k,  such  that 

\ak{v  -  IkPk-iv,v)\  <Ca  a*(u,u)^“",  (3.3.2) 

where  Xk  is  the  maximum  eigenvalue  of 

A. 2:  There  exists  Cr  >  1,  independent  of  k,  such  that 

l|w|P 

— —  <  CR(RkU,u),\/u  e  14, 

A* 

where  Ru  =  {I  —  KlKu)A'^^  and  Ku  =  I  —  RkAu  and  is  an  adjoint  of  with  respect  to  aui-,-)- 

Theorem  3.4.  Let  Bj.  he  defined  by  Algorithm  2.1.  Suppose  that  A.l  and  A. 2  hold.  Let  p  =  2  and 
m(k)  =  TO  for  all  k.  Then,  for  m  suffieiently  large,  there  is  a  eonstant  M  independent  of  m  sueh  that 

\ak({I  -  BkAk)u,u)\  <  SAk(u,u),  Mu  G  14, 


with 

— • 

-  M  +  m°‘ 

Remark  3.5.  Theorem  3.1  shows  that  the  eonvergenee  faetor  of  the  W-eyele  multigrid  algorithm  does 
not  depend  the  number  of  levels  if  the  multigrid  method  is  eonvergent. 

The  following  result  concerns  the  variable  R-cycle  multigrid  algorithm. 

Theorem  3.6.  Let  Bj.  he  defined  by  Algorithm  2.1.  Suppose  that  A.l  and  A. 2  hold.  The  number  of 
smoothings,  m(k),  inereases  as  k  deereases  in  sueh  a  way  that 

(3o'm{k)  <  m(k  —  1)  <  j3im{k) 
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with  1  <  Po  <  Pi-  Let  p  =  1.  Then 


r]oak{u,u)  <  ak{BuAuu,u)  <  r]iak{u,u),  \fu  G  14 

holds  with  rjo  >  M  >  0. 

Remark  3.7.  The  eonvergenee  ofV-eyele  multigrid  algorithm  is  dijjfieult  to  show.  However,  in  many 
numerieal  experiments,  the  V-eyele  multigrid  algorithm  also  eonverges  if  the  number  of  smoothings  m  is 
suffieiently  large. 

4.  Numerical  results.  In  this  section,  we  consider  two  second-order  elliptic  problems  and  compare  the 
convergence  behavior  of  multigrid  algorithms  with  a  covolume-based  intergrid  transfer  operator  and  other 
intergrid  transfer  operators  on  Pi -nonconforming  multigrid  algorithms. 

Example  1.  We  consider  the  Laplace  equation  on  the  unit  square 


where 


—Au  =  /,  in  fl, 
r.  du 

l3u  +  —  =g,  on  dn, 
on 


(4.4.1) 


0  at  y  =  0, 1 
1  at  a;  =  0, 1 


and  the  coarsest  primal  triangulation  of  fl  is  shown  Figure  3(a). 

In  Figure  3,  we  show  some  primal  triangulations  which  generated  by  the  nested-mesh  subdivision  method. 
We  compare  convergence  rates  of  LF-cycle  multigrid  methods  with  with  intergrid  transfer  operator  1^,  inter¬ 
grid  transfer  operator  ij.,  which  was  defined  in  [2],  and  intergrid  transfer  operator  7^^,  which  was  defined  in 
[11]  in  Figure  4(a)  and  4(b).  We  perform  3  Gauss-Seidel  sweeps  in  Figure  4(a)  and  4  damped  Jacobi  sweeps 
with  damping  coefficient  to  =  0.6  in  Figure  4(b).  These  figures  show  that  the  error  reduction  of  is  very 
poor  for  the  Robin  or  the  Neumann  boundary  value  problem. 

We  report  the  average  error  reduction  of  the  LF-cycle  multigrid  method  with  the  covolume  based  intergrid 
transfer  operator  with  respect  to  number  of  grid  levels  in  Table  1(a)  for  the  Gauss-Seidel  smoother  and  Table 
1(b)  for  the  damped  Jacobi  smoother.  For  reference,  we  report  that  the  average  error  reduction  of  the  R-cycle 
multigrid  method  Table  1(c)  and  (d). 

Remark  4.1.  Table  1  shows  that  the  average  error  reduetion  of  the  W-eyele  multigrid  method  does  not 
depend  the  number  of  grid  levels  but  that  of  the  V-eyele  multigrid  method  does  depend  on  the  number  of 
levels.  Also,  Table  1  shows  that  the  V-eyele  multigrid  method  eonverges  for  a  suffieient  number  of  smoothing 
steps. 

In  Figure  5,  we  compare  convergence  rates  of  the  Preconditioned  Conjugate  Gradient  method  with 
multigrid  as  a  preconditioner.  Also,  we  report  in  Table  2  the  average  error  reduction  of  PCG  with  Pi 
nonconforming  multigrid  algorithm  as  a  preconditioner. 

Remark  4.2.  Table  2  shows  that  the  average  error  reduetion  of  PCG  with  variable  V-eyele  multigrid  as 
a  preeonditioner  does  not  depend  the  number  of  grid  levels,  but  a  V-eyele  multigrid  algorithm  does  depend 
on  the  number  of  levels. 

Example  2.  We  consider  the  second-order  partial  differential  equation  with  discontinuous  coefficients 
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Table  1.  Average  error  reduetions  of  multigrid  method,  (a)  W-eyele  multigrid  with  Gauss-Seidel  smoother, 
(b)  W-eyele  multigrid  with  Jaeohi  smoother,  (e)  V-eyele  multigrid  with  Gauss-Seidel  smoother,  (d)  V-eyele 
multigrid  with  Jaeohi  smoother.  (*  :  not  eonverge) 


J\m 

1 

2 

3 

4 

5 

3 

0.667 

0.391 

0.294 

0.185 

0.135 

4 

0.657 

0.386 

0.300 

0.179 

0.137 

5 

0.677 

0.377 

0.300 

0.173 

0.137 

6 

0.677 

0.367 

0.300 

0.165 

0.136 

7 

0.677 

0.360 

0.300 

0.160 

0.136 

(a) 


J\m 

3 

4 

5 

6 

7 

3 

0.420 

0.251 

0.168 

0.152 

0.093 

4 

0.496 

0.264 

0.184 

0.161 

0.105 

5 

0.572 

0.272 

0.196 

0.166 

0.112 

6 

0.646 

0.276 

0.204 

0.169 

0.117 

7 

0.720 

0.281 

0.214 

0.174 

0.121 

(c) 


J\m 

1 

3 

5 

7 

9 

3 

0.778 

0.492 

0.325 

0.221 

0.153 

4 

0.781 

0.489 

0.319 

0.213 

0.146 

5 

0.782 

0.484 

0.311 

0.207 

0.140 

6 

0.782 

0.482 

0.306 

0.201 

0.135 

7 

0.782 

0.480 

0.303 

0.197 

0.131 

(b) 


J\m 

4 

5 

6 

7 

8 

3 

0.535 

0.416 

0.346 

0.291 

0.246 

4 

0.681 

0.503 

0.378 

0.308 

0.261 

5 

0.825 

0.592 

0.432 

0.326 

0.270 

6 

0.973 

0.677 

0.484 

0.353 

0.278 

7 

* 

0.761 

0.535 

0.384 

0.288 

(d) 


T  tI  tII 

Ik —  Ik'  '  h 


(a) 


(b) 

Fig.  5.  Gomparison  of  error  reduetion  of  Preeonditioned  GG.  (a)  Using  1  Gauss-Seidel  smoothing  step  per 
level,  (b)  Using  2  Damped  Jaeohi  smoothing  steps  with  to  =  0.6  per  level 
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Table  2.  Average  error  reduetions  of  Preeonditioned  CG.  (a)  Variable  V-eyele  multigrid  with  Gauss-Seidel 
smoother,  (b)  Variable  V-eyele  multigrid  with  Jaeobi  smoother,  (e)  V-eyele  multigrid  with  Gauss-Seidel 
smoother,  (d)  V-eyele  multigrid  with  Jaeobi  smoother 


J\m{J) 

1 

2 

3 

4 

3 

0.380 

0.175 

0.147 

0.077 

4 

0.387 

0.179 

0.154 

0.078 

5 

0.392 

0.181 

0.156 

0.078 

6 

0.395 

0.181 

0.156 

0.078 

7 

0.397 

0.182 

0.157 

0.077 

(a) 

J\m 

1 

2 

3 

4 

3 

0.490 

0.226 

0.205 

0.104 

4 

0.546 

0.242 

0.229 

0.110 

5 

0.593 

0.256 

0.249 

0.116 

6 

0.634 

0.268 

0.265 

0.121 

7 

0.668 

0.278 

0.280 

0.124 

(c) 


J\m{J) 

1 

2 

3 

4 

3 

0.510 

0.358 

0.262 

0.202 

4 

0.526 

0.362 

0.267 

0.206 

5 

0.535 

0.366 

0.267 

0.207 

6 

0.540 

0.367 

0.269 

0.208 

7 

0.542 

0.368 

0.270 

0.207 

(b) 


J\m 

1 

2 

3 

4 

3 

0.563 

0.433 

0.347 

0.281 

4 

0.615 

0.475 

0.387 

0.312 

5 

0.660 

0.521 

0.420 

0.340 

6 

0.696 

0.559 

0.451 

0.363 

7 

0.732 

0.590 

0.476 

0.385 

(d) 


on  the  unit  square 


—V  •  AVu  =  /,  in 

jJu  +  A—  =  o,  on  dQ., 
on 


(4.4.2) 


where 


A{x,y) 


10,  if  1/3  <  a;  <  2/3  and  1/3  <  y  <  2/3, 
1,  otherwise  , 


P  is  the  same  value  as  in  Example  1  and  the  coarsest  primal  triangulation  of  fl  is  shown  Figure  6(a). 

In  Figure  6,  we  show  some  primal  triangulations  which  generated  by  the  nested-mesh  subdivision  method. 

In  Figure  7  and  8  and  Table  3  and  4,  we  report  the  same  numerical  experiments  applied  to  (4.4.2).  In 
these  numerical  experiments,  the  error  reduction  factor  slightly  increases  in  the  IF-cycle  multigrid  method 
and  PCG  with  variable  F-cycle  multigrid  preconditioner,  but  is  not  rapidly  increasing  compared  with  V- 
cycle  multigrid  and  PCG  with  a  F-cycle  multigrid  preconditioner.  This  result  is  related  to  the  regularity  of 
the  partial  differential  equation. 

Example  3.  As  a  nonlinear  example,  we  consider  a  non-equilibrium  radiation  diffusion  equation  system, 
which  can  be  written  as 


with 


—  -S/-{DrVE)=aa{T^-E),  in  fl  , 
QT 

—  -S/-{DtVT)  =  -aa{T^-E),  in  fl  , 


Oa 


’ 


Dr{T,E) 


1 

3<7„  + A|VE|’ 


DtiT)  =  KTi. 


(4.4.3) 
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Here,  E  represents  the  photon  energy,  T  is  the  material  temperature,  and  k  is  the  material  conductivity. 
In  the  non-equilibrium  case,  the  nonlinear  source  terms  on  the  right-hand-side  are  nonzero  and  govern 
the  transfer  of  energy  between  the  radiation  field  and  material  temperature.  Additional  nonlinearities  are 
generated  by  the  particular  form  of  the  diffusion  coefficients,  which  are  functions  of  the  E  and  T  variables.  In 
particular,  the  energy  diffusion  coefficient,  Dr{T,E)  contains  the  term  |Vi?|  which  refers  to  the  gradient  of 
E.  This  limiter  term  is  an  artificial  means  of  ensuring  physically  meaningful  energy  propagation  speeds  (i.e. 
no  signal  speeds  faster  than  the  speed  of  light)  ([16,  13,  15]).  The  atomic  number  2:  is  a  material  coefficient, 
and  while  it  may  be  highly  variable,  it  is  a  function  of  position  only  (i.e.  2  =  f{x,y)  in  two  dimensions). 

Equations  (4.4.3)  represent  a  system  of  coupled  nonlinear  partial  differential  equations  which  must  be 
discretized  in  space  and  time.  The  time  derivatives  are  discretized  as  first-order  backwards  differences,  with 
lumping  of  the  mass  matrix,  leading  to  an  implicit  scheme  which  requires  the  solution  of  a  nonlinear  problem 
at  each  time  step.  This  approach  is  first-order  accurate  in  time,  and  is  chosen  merely  for  convenience,  since 
the  principal  objective  is  the  study  of  the  solution  of  the  nonlinear  system.  Spatial  discretization  on  two- 
dimensional  triangular  meshes  is  achieved  by  a  Pi-  nonconforming  finite  element  procedure,  assuming  linear 
variations  of  E  and  T  over  a  triangular  element. 

The  test  case  chosen  for  this  work  is  taken  from  [16,  15],  we  consider  a  unit  square  domain  of  two 
dissimilar  materials,  where  the  outer  region  contains  an  atomic  number  of  2  =  1  and  the  inner  region 
(1/3  <  X  <  2/3, 1/3  <  y  <  2/3)  contains  an  atomic  number  of  2  =  10.  The  top  and  bottom  walls  are 
insulated,  and  inlet  outlet  boundaries  are  specified  using  mixed  Robin  boundary  conditions,  as  shown  in  the 
Figure  9.  We  use  the  same  triangulation  to  the  previous  example  2. 
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(a) 


(b) 

Fig.  7.  Comparison  of  error  reduction  of  the  W -cycle  multigrid  method,  (a)  Using  3  Gauss-Seidel 
smoothing  steps  per  level,  (b)  Using  4  Damped  Jacobi  smoothing  steps  with  to  =  0.6  per  level 


Table  3.  Average  error  reductions  of  multigrid  method,  (a)  W -cycle  multigrid  with  Gauss-Seidel  smoother, 
(b)  W -cycle  multigrid  with  Jacobi  smoother,  (c)  V -cycle  multigrid  with  Gauss-Seidel  smoother,  (d)  V -cycle 
multigrid  with  Jacobi  smoother.  (*:  not  converge) 


J\m 

1 

2 

3 

4 

5 

3 

0.783 

0.527 

0.498 

0.311 

0.330 

4 

0.801 

0.549 

0.525 

0.343 

0.357 

5 

0.812 

0.557 

0.539 

0.350 

0.366 

6 

0.818 

0.560 

0.545 

0.355 

0.372 

(a) 


J\m 

3 

5 

7 

9 

11 

3 

0.564 

0.358 

0.245 

0.169 

0.120 

4 

0.722 

0.403 

0.285 

0.203 

0.145 

5 

0.886 

0.473 

0.313 

0.226 

0.164 

6 

* 

0.680 

0.409 

0.264 

0.182 

(c) 


J\m 

1 

3 

5 

7 

9 

3 

0.857 

0.646 

0.511 

0.412 

0.340 

4 

0.872 

0.675 

0.542 

0.445 

0.371 

5 

0.880 

0.689 

0.553 

0.454 

0.377 

6 

0.885 

0.696 

0.557 

0.452 

0.375 

(b) 


J\m 

4 

6 

8 

10 

12 

3 

0.759 

0.492 

0.394 

0.326 

0.270 

4 

0.967 

0.624 

0.447 

0.376 

0.316 

5 

* 

0.751 

0.501 

0.408 

0.346 

6 

* 

0.874 

0.567 

0.432 

0.366 

(d) 
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0  10  20  30  40  50  60 

iterations 


(a) 


0  10  20  30  40 

iterations 


(b) 

Fig.  8.  Comparison  of  error  reduction  of  Preconditioned  CG.  (a)  Using  1  Gauss-Seidel  smoothing  step  per 
level,  (b)  Using  2  Damped  Jacobi  smoothing  steps  with  to  =  0.6  per  level 


Table  4.  Average  error  reductions  of  Preconditioned  CG.  (a)  Variable  V -cycle  multigrid  with  Gauss-Seidel 
smoother,  (b)  Variable  V -cycle  multigrid  with  Jacobi  smoother,  (c)  V -cycle  multigrid  with  Gauss-Seidel 
smoother,  (d)  V -cycle  multigrid  with  Jacobi  smoother 


2  3 

0.299  0.268 

0.311  0.282 

0.314  0.293 

0.315  0.301 

(a) 


J\m 

1 

2 

3 

4 

3 

0.579 

0.327 

0.305 

0.191 

4 

0.632 

0.363 

0.351 

0.219 

5 

0.676 

0.388 

0.391 

0.231 

6 

0.716 

0.412 

0.439 

0.245 

(c) 


4 

0.161 

0.173 

0.178 

0.178 


J\to(  J) 

1 

2 

3 

4 

3 

0.614 

0.465 

0.372 

0.309 

4 

0.647 

0.499 

0.405 

0.343 

5 

0.674 

0.527 

0.429 

0.357 

6 

0.690 

0.541 

0.442 

0.371 

(b) 


J\m 

1 

2 

3 

4 

3 

0.641 

0.500 

0.413 

0.351 

4 

0.699 

0.558 

0.469 

0.408 

5 

0.742 

0.606 

0.515 

0.449 

6 

0.777 

0.643 

0.554 

0.486 

(d) 
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Fig.  9.  Test  problem  for  radiation  transport  problem 


Fig.  10.  Some  numerieal  results  of  radiation  transport  problem 

This  problem  is  highly  nonlinear  and  has  been  identified  as  one  of  the  most  time-consuming  components 
in  large  multiphysics  simulation  codes.  To  solve  this  nonlinear  problem,  we  use  Newton  linearization  method. 
However  the  resulting  linear  problems  are  nonsymmetric,  so  we  use  Preconditioned  GMRES(PGMRES)  with 
a  nonconforming  multigrid  preconditioner  (F-cycle  or  the  variable  F-cycle  multigrid  algorithm  defined  in 
section  3)  to  solve  the  linear  problem.  In  Figure  10,  we  illustrate  a  typical  simulation  result  for  this  system. 
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We  plot  the  contour  of  temperature  T  at  time  t  =  2.0,  3. 0.4.0,  5.0.  These  show  that  the  solutions  are  rapidly 
changing  near  the  position  where  the  two  different  materials  meet. 

To  get  the  coarse  grid  operators,  we  need  to  get  some  approximations  in  the  coarse  finite  element  space 
of  solutions  in  the  finest  finite  element  space.  When  we  use  {iD^  as  the  fine-to-coarse  intergrid  transfer 
operator,  we  cannot  solve  some  coarse  level  problems  because  coarse  level  bilinear  operator  fails  to  be 
defined  (some  values  of  temperature  T  are  negative,  as  remark  3.2),  or  the  coarse  level  problem  is  very  hard 
to  solve  by  using  iterative  methods.  But,  if  we  use  (Ik)^  as  fine-to-coarse  intergrid  transfer  operator,  the 
preconditioners  work  well. 

In  Figure  11,  we  compare  the  error  reduction  of  PGMRES  with  E-cycle  and  variable  E-cycle  with 
smoothing  number  1  and  2  and  Gauss-Seidel  smoothing  at  time  t  =  2.0  with  time  step  size  dt  =  0.001,  0.002, 
0.005,  0.01.  We  measure  linear  residual  error  by  preconditioned  error  and  stop  the  linear  PGMRES  iteration 
if  the  relative  linear  residual  error  is  less  than  10“®.  Eor  each  time  step,  we  stop  the  nonlinear  iteration  if 
the  nonlinear  residual  error  is  less  than  5  x  10“®. 

The  numerical  results  at  time  t  =  2.0  show  that  there  is  no  significant  difference  between  E-cycle  and 
variable  E-cycle  preconditioner,  but  significant  improvement  in  error  reduction  in  solving  the  linear  problem 
when  one  increases  the  smoothing  number. 

We  have  a  studied  covolume-based  intergrid  transfer  operators  in  a  Pi  nonconforming  multigrid.  We 
found  that  multigrid  methods  with  covolume-based  intergrid  transfer  operators  converge  more  slowly  than 
with  a  standard  previous  intergrid  transfer  operator.  This  result  was  expected  because  this  operator  is  simple 
and  does  not  preserve  as  many  high-order  functions  as  the  standard  operator.  However,  we  promote  this 
operator  for  preservation  of  positivity  in  solving  nonlinear  problems  and  for  parallelization. 

Acknowledgments.  The  author  would  like  to  thank  D.  E.  Keyes  of  Old  Dominion  University  for  his 
valuable  advice  in  the  preperation  of  this  paper. 
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